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Abstract: In this paper we discuss the use of wavelet bases to solve the relativistic three-body 
problem. Wavelet bases can be used to transform momentum-space scattering integral equations into 
an approximate system of linear equations with a sparse matrix. This has the potential to reduce 
the size of realistic three-body calculations with minimal loss of accuracy. The wavelet method leads 
to a clean, interaction independent treatment of the scattering singularities which does not require 
any subtractions. 
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I. INTRODUCTION 

This is the third paper 00 in a series of investiga- 
tions designed to explore potential advantages of using 
wavelet numerical analysis to solve the relativistic three- 
body problem. Commercially, wavelets are used to con- 
vert raw digitized photographic images to compressed 
JPEG files 1^]. In this application the data compression 
leads to a large saving in storage space with a minimal 
loss of information. The compression involves expanding 
the raw digital image in a wavelet basis and setting the 
smaller expansion coefficients to zero. The kernel of a 
scattering integral equation and a raw digital image can 
both be approximated by rectangular arrays of numbers 
with some continuity properties. This suggests that the 
bases used to compress digital images could be used to 
generate accurate sparse matrix approximations to the 
kernel. 

The ability to construct numerically exact solutions 
to the quantum mechanical three-body problem cou- 
pled with the ability to accurately measure complete 
sets of experimental observables constrains the form of 
the three-nucleon Hamiltonian. These constraints have 
resulted in the construction of realistic model nucleon- 
nucleon interactions 0[3 0- When these interactions 
are used in the many-nucleon Hamiltonian, the resulting 
dynamical model provides a good quantitative descrip- 
tion of low-energy nuclear physics Q • 

The state of the art in few-body computations has im- 
proved to the point where numerically exact scattering 
calculations at energy and momentum transfers of hun- 
dreds of MeV have been performed Q. Higher energy 
calculations are possible. As in the low-energy case, the 
structure of Hamiltonians for higher energy reactions can 
be constrained by the consistency of the few-body cal- 
culations with precise measurements of complete sets of 
experimental observables. 
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The success of the few-body approach to low-energy 
nuclear physics is a consequence of (1) knowing the rel- 
evant degrees of freedom (nucleons), (2) working with 
the most general Hamiltonians involving these degrees 
of freedom that are consistent with the symmetries of 
the system (Galilean invariance) and (3) understanding 
the relation between the few and many-body problem 
(cluster properties). To extend this success to reactions 
involving higher energy scales (1) the relevant degrees 
of freedom may have to include explicit mesonic or sub- 
nucleonic degrees of freedom (2) Galilean invariance must 
be replaced by Poincare invariance and (3) cluster prop- 
erties must be maintained. 

Each of the required extensions of low-energy nuclear 
dynamics is non-trivial, and progress has been made on 
all three problems JJ jE],10],ll] .12j- The purpose of this 
paper is to focus on technical aspects of using wavelet 
numerical analysis to construct exact numerical solutions 
of the dynamics in Poincare invariant few-body models. 
While the scope of this paper is limited to three-nucleon 
models with no explicit mesonic or subnucleonic degrees 
of freedom and S- matrix cluster properties p^ . the for- 
mulation and advantage of methods discussed in this pa- 
per are straightforward to extend to systems with explicit 
mesonic degrees of freedom and stronger forms of cluster 
properties Il2l| . Moving singularities are a generic feature 
of the dynamical equations in all of these cases. 

Relativistic few-body equations are naturally formu- 
lated in momentum space. Relativistic kinematic factors, 
Wigner rotations and Melosh rotations are all multiplica- 
tion operators in momentum space. The compactness of 
the iterated Faddeev-Lovelace kernel implies the kernel of 
the integral equations can be uniformly approximated by 
a finite matrix, resulting in a finite linear system. In the 
momentum representation these linear systems have large 
dense matrices, which increase in size with increasing en- 
ergy and momentum transfer. It is desirable to be able to 
perform accurate calculations at energy and momentum 
scales where subnuclear degrees of freedom are relevant. 
At these scales a relativistic treatment is required and ad- 
vances in computational efhciency are needed to perform 
realistic calculations. The ability of the wavelet trans- 
form to efficiently transform a dense matrix to an ap- 
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proximate sparse matrix suggests that wavelet methods 
can provide a powerful tool for improving the efhciency 
of relativistic few-body computations. 

The advantages of using wavelet numerical analysis to 
solve momentum space scattering integral equations was 
investigated in Q and 0. These papers used wavelet nu- 
merical analysis to solve the Lippmann-Schwinger equa- 
tion for a system of two nucleons interacting with a 
Malfliet Tjon V potential 0| using partial wave ex- 
pansions 1] and direct integration In both applica- 
tions the kernel of the integral equation was accurately 
approximated by a sparse matrix, which resulted in ac- 
curate approximate solutions. The success of these ap- 
plications indicates that wavelet numerical analysis will 
have similar advantages when applied to the relativistic 
three-body problem. 

The feature of the three-body problem that is not 
present in the two-body applications is moving singular- 
ities. The methods used in ^ and '3| are not applicable 
to problems with moving singularities. The purpose of 
this paper is to illustrate how to apply wavelet numerical 
analysis to treat the moving singularities that appear in 
the relativistic three-body problem. 



II. OVERVIEW - WAVELET NUMERICAL 
ANALYSIS 

The applications in ref. and used Daubechies' 
wavelets. These wavelets differ from the wavelets used 
to store JPEG images. The Daubechies' wavelets are 
orthonormal but the basis functions are not reflection 
symmetric; while the wavelets used to store JPEG im- 
ages sacrifice orthonormality to obtain more symmetric 
basis functions. The results of ref. |l| and [J indicate 
that the Daubechies' wavelets are suitable for scattering 
calculations. 

Daubechies' wavelets are discussed in many 

texts on wavelets jl^j JJJ ^20|| They are fractal 

functions that have complex structures on all scales. Be- 
cause the basis functions have structure on all scales, 
numerical applications with wavelets require a different 
approach to numerical analysis, hence the term wavelet 
numerical analysis. 

We use the Daubechies' wavelets because they are a 
dense orthonormal set of compactly supported functions 
with the property that finite linear combinations can lo- 
cally pointwise represent low-degree polynomials. 

Two types of functions are needed to generate wavelet 
bases. These functions are called scaling functions and 
wavelets. The scaling function, <j>{x), is the solution of 
the linear renormalization group equation: 

2K-1 

Dcj^ix) = J2 hiT'H^) (1) 



Table 1: Daubechies' K — 3 Scaling Coefficients 

~ K=3 
Ik) (1 + VW + \/ 5 + 27T0) /16\/2 
hi (5 + VW + 3\/ 5 + 2^/I0)/ 16^/2 
h2 (10 - 2710 + 2\ /5 + 2^ )/16^ 
hs (10 - 2VT0 - 2 ^/5 + 2yi 0)/16^ 
hi (5 + VW - 3\ /5 + 2v/l0 )/16\/2 
hs (1-H VTO- \/5 + 27T0)/16V2 



with normalization 

/oo 
(j>{x)dx = 1. (2) 
-00 

Equation is called the scaling equation. 

In equation Z? is the unitary scaling operator 

Df{x) -^/(|) (3) 

which stretches the support of the function by a factor 
of two. The operator T is the unitary unit translation 
operator 

Tf{x) = f{x - 1). (4) 

The coefficients hi are real numbers that determine the 
properties of the scaling function. X is a finite positive 
integer. The calculations in l] 2] used Daubechies' K — 
3 wavelets. The reason for this choice will be discussed 
later. For the K — 3 Daubechies' wavelets the six scaling 
coefficients hi are given in Table 1. 



The fractal structure of (j){x) is a consequence of the 
scaling equation |^ which shows that the scaling function 
on a given scale is a finite linear combination of translates 
of the same function on half the scale. 

The scaling equation implies that the scaling coeffi- 
cients hi satisfy 

2K-1 

Y,hi^V2 (5) 

and that the solution (t){x) of equation ^ has support 
on the interval [0, 2K - I] p^ . 

The unit translates of the scaling function are or- 
thonormal 

{T'-<j>,T"cf,)^S,r^n (6) 

provided the scaling coefficients satisfy the additional 
constraints: 

2K-1 

hlhi^2m = ^mO- (7) 

Z=0 
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FIG. 1: Daubechies' K = 3 scaling function. 

The scaling function (j){x) is continuous (for K > 1) 
and can be computed exactly at all dyadic rationals using 
equations ^ and Q. This method is used to compute 
the Daubechies' K — 3 scaling function plotted in Figure 
1. 

The subspace of square integrable functions on the real 
line that can be expressed as linear combinations integer 
translates of the scaling function, T"0(x), is the subspace 
Vo of L^{R) defined by : 

OC OO 

Vo:={/(x)= /"^^"'^(^)l E l/"l'<°«}- (8) 

n— — oo n— — oo 

Application of powers of the scaling operator D'^ to 
Vo defines subspaces Vk with coarser {k > 0) or finer 
resolution {k <0): 

Vk - D'^Vo. (9) 

The space Vk is called the approximation space with reso- 
lution k. The resolution determines the size of the small- 
est features that can be approximated by functions in 
Vk- 

The scaling functions 

q^Ux) ■■= D^r^m - J - n) (10) 

are an orthonormal basis for Vk- The support of (t>kn{x) 
is [2''n,2''{n + 2K 

The scaling equation implies the inclusions 

Vk^Vk+i. (11) 

The orthogonal compliment of Vfe+i in Vk is denoted by 
Wk+i which leads to the orthogonal decomposition 

Vk = Vk+i®Wk+i. (12) 




X 



FIG. 2: Daubechies' K — 3 mother wavelet. 

Orthonormal basis functions, ^pkm{x), for the sub- 
spaces yVfc are elements of Vk-i given by 

^km{x) = D'^r^i^ix) (13) 



2K-1 

^(x) = 9iD''T'^{x) (14) 

1=0 

where 

gi = {-yh2K^i^i. (15) 

The subspaces Wkn are called the wavelet spaces and the 
basis functions -ipkn are called wavelets. The support of 
ipknix) is identical to the support of 4>kn{x)- Since the 
wavelets are finite linear combinations scaling functions 
they are also fractal functions. 

The function il){x) = ipooix) is called the mother 
wavelet; the Daubechies' K = 3 mother wavelet is shown 
in Figure 2. 

The coefficients hi for the Daubechies' K wavelets are 
determined by equations lO and {Tj) and the require- 
ments 

/OO 
^p{x)x''dx = 0; n = Q,l,--- ,K-1 (16) 
-OO 

which implies that ipkmi^) is locally orthogonal to all 
polynomials of degree K — 1. Equation (|16|) can be ex- 
pressed directly in terms of the scaling coefBcients: 

2K-1 
1=0 
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2K-1 



l\-yh2K-i-i^0 k^O---K-l. (17) 



For any K > Q conditions |(SJ|, 10 and l(T7)l determine hi 
up to reflection, 



hi — > h'l :— h2K-i- 



(18) 



The entries in Table 1 are the solution of these equations 
for = 3. The Daubechies' wavelets have the property 
that as fc — > — CX3 (infinitely fine resolution) the space Vk 
becomes all of L^(R). 

The decomposition H12(l implies that 

Vfc = Wfc+i e Wfe+2 • • • Wfe+m-i e Wk+,n © Vfe+™, (i9) 

for any m > 0. This means that functions in the ap- 
proximation space Vk can be expanded as linear combi- 
nations of the scaling basis functions of resolution k or 
equivalently as linear combinations of the scaling basis 



functions of a coarser resolution k' = k 



and wavelet 



basis functions of all resolutions from fc + 1 to fc + m. 

If we let A: — > — cxD with m + k finite, then the functions 
in H19|) become a basis for L^(IR). Since the i/jknix) basis 
functions are locally orthogonal to degree K — 1 polyno- 
mials, and only a finite number of the (j)kn (x) are non- 
zero at any x, it follows that finite linear combinations of 
4>k,n{x) must be able to locally pointwise represent degree 
K ~ \ polynomials. 

Thus, for the Daubechies' K — 3-wavelets finite lin- 
ear combinations of the scaling basis functions (t>kn{x) 
can locally pointwise represent polynomials of degree 2, 
while the wavelet basis functions tpkn{x) are orthogonal 
to degree 2 polynomials. 

Equation 1)19(1 implies that the projection Pk of a func- 
tion f{x) on Vk can be represented by 



Pfc/(x) = ^a„(/)fe„(a;) a„ = j f(x)(j3kn{x)dx (20) 

n •' 

or equivalently 

PkJjx) = y^a„(/)fc+r«,«(a^) + bk'n'ipk'^n, (21) 



k'=k+l 



bkn = / f{x)ll}k,n{x)dx. 



(22) 



For a sufficiently fine resolution (large — fc) the scal- 
ing basis functions have small support and integrate to 
a constant. If f{x) varies slowly on intervals of width 
(2K — 1)2^ then expansion coefficients a„ are well ap- 
proximated by evaluating f{x) at any point in the sup- 
port of (j)kn (x) and multiplying by 2^"'/^ , which is the inte- 
gral of (j)kn {x) ■ This means that the scaling function ba- 
sis coefficients, a„, are well approximated, up to a fixed 



multiplicative constant, by sampling the original func- 
tion. These coefficients play the role of the raw image 
in a digital photograph. They provide an accurate, but 
inefficient approximation of the function f{x). 

In the representation (|21|l . if f(x) can be accurately 
approximated by a polynomial of degree K — 1 on the 
support of ipknix) then bkn ~ 0. This means that if f{x) 
can be well approximated by low-degree local polynomi- 
als with compact support on multiple scales, most of the 
coefficients bkn will be be small and the function can be 
accurately approximated by replacing these small coeffi- 
cients by zero. The mean square error is the sum of the 
squares of the discarded coefficients, which can be con- 
trolled by selecting a maximum size of the discarded co- 
efficients. Even though most of the basis functions in the 
representation H21|) are orthogonal to low degree polyno- 
mials, the equivalence between the representations (|20() 
and (EU means that the representation l(21() can still lo- 
cally pointwise represent low degree polynomials. The 
orthogonal transformation connecting equivalent repre- 
sentations (|20|l and (|21|l of Vk is called the wavelet trans- 
form |2||. For N basis functions it can be implemented 
in 0{N) steps, which for large TV requires less steps than 
a fast Fourier transform. 

The scaling equation and normalization condition can 
be used to derive exact expressions for the moments and 
partial moments of the scaling function 



(j)kn{x)x"'dx 



(23) 



(j)kn ix)x''' 



< LI' <2K -1 



2''l 



(24) 

in terms of the scaling coefficients hi . Explicit expressions 
for the moments and partial moments appear in 0|, 

M 

For the Daubechies' K wavelets with K > 1 the second 
moment of the scaling function is the square of the first 
moment. This means that for K = Z the first moment 
provides a single quadrature point that will integrate the 
scaling function times any second degree polynomial ex- 
actly: 



(j){x){a -I- 6a; -I- cx^)dx — a + b{x) + c{x)'' 



(25) 



This is called the one-point quadrature. Translating and 
rescaling leads to one-point quadratures for all of the 
scaling basis functions (j>kn{x). The choice K — 3 
Daubechies' wavelets in [J] and 0| is motivated by their 
ability to locally pointwise represent second degree poly- 
nomials and to exactly integrate these local polynomials 
with a one-point quadrature. 

In 0, and the moments and the scaling equa- 
tion are used to compute the singular integrals 



kn 



— dx 



(t'kn jx) 

X zLie 



(26) 
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to any pre-determined precision. 

In applications the integral equation is approximated 
by projecting on an approximation space Vk with a finest 
resolution k dictated by the problem. This projection can 
be computed efficiently in the scaling basis (|20|l using the 
one-point quadrature H25|l and the explicit integrals (|26|l . 
The resulting matrix equation is transformed using the 
wavelet transform to an equivalent system in the wavelet 
basis (|21|) . In the transformed basis the kernel of the 
integral equation decomposes into the sum of a sparse 
matrix and a small matrix. The kernel is approximated 
by setting matrix elements of the kernel that are smaller 
than a threshold value to zero. This results in a sparse 
matrix approximation. The resulting linear system is 
solved using sparse matrix iterative techniques, such as 
the complex biconjugate gradient method 23i||2J] used in 
. This solution is transformed back to the scaling func- 
tion representation, using the inverse wavelet transform, 
and the resulting solution is inserted back in the integral 
equation to construct an interpolated solution psj . 

The only wavelet information used in this application 
is the wavelet transform and moments of the scaling func- 
tion. These can both be expressed directly in terms of 
the scaling coefficients hi in Table 1. This means that 
the basis functions never have to be computed. 

The work in references [2| shows that all of these 
steps work as expected. These references also discuss 
technical issues that arise due to the treatment of end- 
points when the equations are transformed to a finite 
interval. 



III. DYNAMICAL EQUATIONS 

The general structure of the Faddeev-Lovelace [T^|^ 
equation in a relativistic quantum theory with three par- 
ticles of mass m is 



X{p,q;p',q')=D{p,q;p',q'y 



K{p,q-p",q",z)dp"dq" 
z-ei{p",q")-e2{q") 



Xip",q";p',q') (27) 



where 



ei(jj,q) — \/ 4p2 -|- 47n2 ^ ^2 62(5) = \/ rn? -t- q^ 

(28) 

and X ~ X„in, K — Kmn and D = Z3„m are complex 
matrix valued functions. The quantity K(j), q]p' , q' , z) is 
the smooth part of the kernel. 

In order to use wavelet methods it is advantageous to 
transform this equation to a form where functions of the 
momentum, rather than the energy, are additive in the 
denominator. This transformation facilitates the treat- 
ment of the moving singularity. Note that ei > 62 for all 
values of p and q. If i? > then E + ei{p, q) — e2{q) > 0. 
It follows that the singular denominator 

i? + zO+-ei(p,g)-e2(g)' ^^^^ 



where z = E + iO"*", can be transformed to a more useful 
form by multiplying the numerator and denominator by 
the non-zero function E + ei(p, q) — 62(9). This leads to 
the equivalent expression 

1 



E + iO+ -ei- 62 



E + ei — 62 



-h e| - 2Ee2 - ej + iO+{E + ei - €2) 



E + ei — 62 



E^ - 3m2 - 2Ey/q^ + - V + iO+ 



(30) 



which has the advantage that it separates the p and q de- 
pendence. In this expression there is only one singularity 
in the denominator. The next step is to change variables 



X = rjAp"^ y — r]lE(^\J (p- + vn? — m) (31) 
and define 



z' = Tj[{E ~ mf - 4m^]. 



(32) 



The parameter rj both sets a scale and can be used to fine 
tune z' so the real part is a dyadic rational of the form 
n/2~''. The method that we use to evaluate the singular 
integrals requires that z is a dyadic rational. 

The substitutions H31|l and H32|) lead to the equivalent 
equation 

X{x, y; x', y') = D{x, y; x' , y') + 



Kix,y;x'\y")dx"dy" 
z' - x" - y" + iO+ 



X{x",y";x',y') (33) 



where 



Xix, y; x', y') = X{p{x),q{y)-p{x'),q{y')) (34) 



D{x, y: x', y') = D{p{x), q{y):p{x'),q{y')) (35) 



Kix,y;x ,y) ^ K{p{x),q{y);p{x ),q{y ),z)x 



^[E + 6,{p{x'),q{y'))-62{q{y')))\^^,\. (36) 

Approximate equations are derived using projection 
methods. We seek a solution X, in the x, and y vari- 
ables, on the approximation space Vfe x V^. Approximate 
equations are obtained by projecting the smooth part of 
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the kernel and the driving on this space. We use the 
approximations: 

^{x,y;x',y') w 



'^4>kTn{x)(l)kn{y)Xm,n{x' ,y'), (37) 



and 



D{x,y\x',y') 



^ (f>hm{x)(j>kn{y)Drn,n{x' ^y')^ (38) 



K{x,y]x',y' 



4'km{x)(t)kn{y)Km {x')4>kn'{y') (39) 

mnm' 7i' 

where 



-D™,n(x', ?/') := 2'=7^(a;„, x„, x', y') (40) 



Table 2: Test of double expansion 

l.OOOOOOe + 04 
6.561000e + 03 
4.096000e + 03 
2.401000e + 03 
1.296000e + 03 
6.250002e + 02 
2.560002e + 02 
8.100015e + 01 
l.eOOOlOe + 01 
5.062574e + 00 
1.000050e + 00 
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1 HQ 
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1 nn 
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1 HQ 

+ Uo 
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1 nn 
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i.zyDUUUe 


1 HQ 

+ Uo 


— o.uuuuuue 


1 nn 
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D.zoUUUUe 


1 no 
+ Uz 


— 4.uuuuuue 


1 nn 
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z.oDUUUUe 


1 no 
+ Uz 


— o.uuuuuue 


1 nn 

+ uu 


o.iUUUUUe 


+ Ui 


— z.uuuuuue 


1 nn 
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i.oUUUUUe 


1 ni 
+ Ui 


— i.oUUUUUe 


1 nn 
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o.UDzoUUe 


1 nn 
+ UU 


— i.UUUUUUe 


1 nn 
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i.UUUUUUe 


1 nn 
+ UU 


—o.uuuuuue 


m 
— Ui 


a ocnnnn^i^ 
D.zoUUUUe 


no 
— Uz 


— o. /OUUUUe 


— Ui 


i.y { {Do\je 


no 
— Uz 


— z.ouuuuue 


m 

— Ui 


o.yUDzoUe 


nQ 

— Uo 


— i.zoUUUUe 


— Ui 


z.44i4UDe 


— U4 


u.uuuuuue 


1 nn 
+ UU 


U.uuuuuue 


1 nn 
+ UU 


i.zouuuue 


m 

— Ui 


z.44i4UDe 


f\A 
— U4 


z.ouuuuue 


— Ui 


o.yuDzoue 


nQ 

— Uo 


o. / ouuuue 


m 
— Ui 
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no 
— Uz 


O.uuuuuue 


m 
— Ui 


a oc:nnnn^ 
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i.UUUUUUe 
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i.UUUUUUe 
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+ UU 


z.uuuuuue 
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+ UU 


i.oUUUUUe 


1 ni 
+ Ui 


d.uuuuuue 


+ 00 


O.iUUUUUe 


+ 01 


4.000000e 


+ 00 


2.560000e 


+ 02 


S.OOOOOOe 


+ 00 


6.250000e 


+ 02 


e.OOOOOOe 


+ 00 


1.296000e 


+ 03 


7.000000e 


+ 00 


2.401000e 


+ 03 


S.OOOOOOe 


+ 00 


4.096000e 


+ 03 


g.OOOOOOe 


+ 00 


6.561000e 


+ 03 


l.OOOOOOe 


+ 01 


l.OOOOOOe 


+ 04 



6.252593e 
1.979528e 
3.920094e 
2.519414e 
1.757732e 
2.398553e 
3.895922e 
1.975902e 
6.247759e 
9.999534e 
1.599991e + 01 
8.099986e + 01 
2.559998e + 02 
6.249998e + 02 
1.296000e + 03 
2.401000e + 03 
4.096000e + 03 
6.561000e + 03 
l.OOOOOOe + 04 



02 
02 
03 
04 
06 
04 
03 
02 
02 
01 



^m,n;m' ,n' • ^ ^ \Xrn^ Xnj X^-a' ^ Xn' ) I41j 

are evaluated at the one point quadrature points associ- 
ated with 4>km(x): 

^ 2K-1 

Xjn = 2''(< X > +m) < a; >= ^ Ihi. (42) 

It is useful to approximate the product of the approx- 
imate expressions, (IS7|l and for K{x,y; x' ,y') and 
X{x',y') by expanding the x' and y' dependence in the 
basis on x Vfc. The justification for this approximation 
is that if both of the approximations are well represented 
by low-degree polynomials on the scale k in x' and y' , 
then the product of these functions should be well rep- 
resented by low-degree polynomials on the scale k in x' 
and y'. 

To test this approximation we approximate x^ by 
re-expanding the product of expansions of x'^ using 
Daubechies' K — 3 wavelets. We write 

X^ ^J^^'nhnix) (43) 
n 

which gives 

X'^ = '^xl^xl(t>km{x)(f)kn{x). (44) 
mn 



This is exact for the Daubechies' K — 3 wavelets. We ap- 
proximate this by projecting on the approximation space 
Vfc. The expansion coefficients are 

ci= x^(f)ki{x)dx = ^xl,xll^^i (45) 

mn 

where 

^Imn-^ J <l)kl{x)(j)krnix)(l)kn{x)dx. (46) 

The approximation defined by this projection can be 
written as 

^*«$^a:^4^mn/0M(x) (47) 

mnl 

The expansion coefficients a;^ are computed using the 

1 point quadrature, which is exact for the expansion of 
x^. The (j)ki(x) are evaluated at dyadic rationals so there 
is no error in computing the scaling basis functions. The 
only source of error is the approximation (|47|l . Table 

2 compares the right and left sides of equation (|47|) for 
resolution fc = — 5 
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The expansion is essentially exact, except near the crit- 
ical point, — 0, where it is still accurate. The ac- 
curacy near the critical point can be improved using a 
higher resolution, however a degenerate critical point is 
not generic. 

This additional approximation gives 

Kix,y;x\y')Xix',y';x",y")^ 



)0fc„"'(y')- (48) 

Even though this introduces two additional sums, most of 
the terms are zero because ^, ^„ = unless \m — m'\, 
\m' — m"\ and |m" — m\ are all less than 2K — 1. In 
section IV. we show that the integrals Imm'm" '^^^ tie 
computed analytically using the scaling equation. 

With these approximations the dynamical equations 
reduce to the algebraic system 

XraA^' ,y') = D,n,n{x' ,y') + 



^ ^ -^m,7i;m' ,n' Im"'m'm" -^n"'7i'n" J-m'" ,n"' {^)^m" ,n" , y ) 



where 



(49) 



JU^)--- rdxdy^^t!^^. (50) 



z — X — y + iO+ 



These equations separate the smooth part of the 
physics input in D and K from the singular part 
of this equation, contained in the integrals J^^{z). 
While the construction of the smooth kernel in the 
relativistic case is considerably more complicated than 
in the non-relativistic case |9||l3||. given the driving 
term and smooth kernel the projections Km,n:m' ,n' and 
Dm,n{x' ,y') can be calculated by evaluating the exact 
driving term and kernel at the one-point quadrature 
point for each (j3km{x). This reduces a Galerkin projec- 
tion to a simple function evaluation. 

In the next section we discuss the evaluation of the 
integrals 



[,m,n 



and 



(51) 



that appear in H49|) . These integrals can be evaluated and 
stored before calculation. They are the wavelet input to 
the calculation. They replace all of the integrations in the 
integral equations and they are independent of the choice 
of dynamical model. The physics input is contained in 
the matrices Km,n;m',n' and Dmn{x,y)- Equation H49() 
gives a clean and stable separation of the physics and the 
treatment of the moving singularity, which is contained 
in the integrals H51|) . 



The equations (|49|l are an infinite set of equations. 
They can be reduced to a finite set by including high- 
momentum cutoffs or transforming to a finite inter- 
val. The treatment of endpoints in the evaluations of 
Km,n;m' jj' and Dm,nix' ,y') IS identical to the treatment 
used in |l| and 0, where partial moments of the scaling 
function are used to construct simple quadratures that 
exactly integrate the product of the scaling function and 
degree K — 1 polynomials over a subinterval of the sup- 
port of the scaling function [2^ [l^ ■ The treatment 
of endpoints in the evaluation of and Jm,n{z) is dis- 
cussed in this paper. 

Even with the reduction to a finite set of equations, 
the system of equations is large. It can be reduced by 
performing a wavelet transform on the scaling function 
basis. This can be done following the method used in 
Qll^, which maps the interval to a circle to treat end- 
points. This does not change the final result because the 
resulting transformation is still a finite orthogonal trans- 
formation. 

The next step is to discard the small matrix elements 
in the transformed kernel and to solve the resulting equa- 
tion. 

As discovered in [2], the treatment of the endpoints 
leads to an ill-conditioned matrix. This is because the 
right tail of the scaling function is small (see Fig. 1). 
Some of the overlap integrals with support containing 
the left endpoint replace the orthogonality integrals by 
integrals of products of scaling functions over an interval 
where the product is small. This can be fixed using the 
conditioning method that was used in 2]. The resulting 
conditioned equations are stable and can be accurately 
solved using sparse matrix techniques. 

The resulting solution can be transformed back to the 
scaling function basis. An interpolated solution is then 
constructed from the solution, Xm,n{x,y), of the alge- 
braic equations using the Sloan interpolation method |25j 

X{x, y; x', y') = D{x, y; x' , y')+ 



Km^njx, y)I^^,^,,I^^ij^iiJ^ii ,^ii{z)Xm\n'{x ,y ). 

(52) 

If this interpolation is used the basis functions never have 
to be evaluated. 

Equations H49|) and (|52|) along with the methods for 
computing the integrals H51|l are the main results of this 
paper. 



IV. EVALUATION OF INTEGRALS 



In this section we discuss the evaluation of the integrals 
ij'nm Jmni^) that appear in equation ipi^l . These 
integrals are defined in equations (|46|l and H5U|) . 

To evaluate these integrals we first express the scale 
"A:" integrals in terms of the scale "0" integrals, then we 



8 



evaluate the scale "0" integrals. Using the definition l|10() 
in equations H49I) and H5U|I we obtain 



and 



rk _ n-k/2 tO 

Ln.ni Lvi 



(53) 



(54) 



For k a negative integer, we can choose 2~'^z as an integer 
which is equivalent to choosing z to be a dyadic rational. 
This can be done for any on shell energy by adjusting the 
parameter rj in H31|l . As a result, it is enough to evaluate 
Jm.ni^) ^^'^ ^i m n ^ ™' " integers. In what follows we 
define 



Lm,n 



and 



(55) 



(56) 



Both Ii_n,m and Jm,n{k) involve integrals over the half 
infinite interval. In order to evaluate these integrals we 
first evaluate the corresponding integrals over the infinite 
interval: 



(x — l)(j){x — n)(j){x — m)dx (57) 



and 

Jrn,n{k) = [ dx [ dy^ — in)(j>(yT^^ ^^^^ 



— oo J —oc 



k — x — y + iO+ 



These integrals are easier to compute because of the sim- 
plified boundary conditions. 

The support of the scaling functions implies that if any 
of ^, m or n are non-negative then 

and if any of l,m or n are less than —2K + 2 then 

Il,m,n - 0. (60) 

The non-trivial values of Ii.m,n correspond to the case 
that the indices /, m, n satisfy 



,m,ne [-2K + 2,-1] 



(61) 



To compute the integrals Ii^m,n defined in (|57|) note that 
the definition implies 

Il,m,n = Io,m-l,n-l- (62) 

which allows us to express /; „ „ in terms of /„ „ defined 
by 



dx(/){x)(j){x ~ m)(l>{x ~ n). (63) 



Since the support of is contained in the interval 

[0,2K — 1], there are only a finite number of non-zero 
values of /,„,„. These have to, n e [—2K + 2,2K — 2]. For 
K — 3 there are 81 non-zero Im.n with m,n ^ [—4,4]. 

We can derive linear equations relating these integrals 
using the scaling equation in the form 



2K-1 



0(x) = 72 ^ /ii0(2x - 0- 



(64) 



1=0 



When we use (|64f) in equation (|63|l . the resulting scal- 
ing equations for the integrals Imn are: 



2K-1 

= V2 ^ hi^hi^hi^l2m+l^-lk,2n+l,^-lk- 



(65) 



These are homogeneous equations relating the non-zero 
values of Im.n- An additional inhomogeneous equation 
is needed to solve for the non-zero values of Im.n- The 
needed equation follows from the normalization condition 
(O and the identity 



^(/)(x - n) = 1 



(66) 



which when used in (|63|1 gives the inhomogeneous equa- 
tion 



2K-2 

E ^ 

m=-2K+2 



(67) 



Equations H65|) and (|67|l are a finite system of [AK — 3) x 
[AK — 3) linear equations that can be solved for the non- 
zero values of Imn ■ The results of these calculations Imn 
for K = i are given in Table 3. 

These solutions give when l,m or n are non- 

negative from equations (|59|) and H62|) . To calculate re- 
maining non-zero values of Iimn first observe that using 
(jMjl in gives scaling equations for Ik,m,n' 



V2j2h 



lk^lm^l,zl2k+lk.2m+l,n,2n+l„ 



(68) 



These equations are not homogeneous equations because 
when any of the indices on the right hand side of the 
equation are non-negative, Ik,m,n = Ik,m,n = Im-k,n-k, 
which is known input. This linear system can be solved 
for the non trivial values of Ik.m,n associated with the 
values of k,n,m S [—2K + 2,-1]. For K = 3 there 
are 64 values of k,m,n E [—4,-1]. The results of this 
calculation for the K = 3 case are given in Table 4. 

All of the overlap integrals I^mn ^^at appear in H49|) 
can be computed from the values in the tables using the 
relations (|53|) . H59|) and (|60|l . There are only a finite 
number of these integrals that are non-zero, so they can 
be computed once and stored. 

The second integral that is needed as input to equa- 
tion is Jmn{l)- The first step to compute Jmn{l) 
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Table 3 - /„ 



Table 4 - /„ 
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1 










-3 





-8.249248e 


- 04 


2 


4 


-2.811543e 


-06 


Jii ~ 


■E 


/ i 


' dy- 


tC 




~ yY(j){x)(j){y) 




-2 





2.283264e 


- 02 


3 


4 


6.184412e 


-06 




n 




1 J — OO ^ 


— OQ 


ft" 










-1 







-8.660587e 
9.104482e 


- 02 

- 01 


4 


4 


-4.467813e 


-06 


conver, 


ges 


. Using the binomial theorem we can exp 



is to compute Jmn{l) defined in (jSHl- With a change of 
variables this integral can be rewritten as 



the integrals in this series in terms of the known moments 
of the scahng function 



dx' 



dy' 



Jk- 



H^'My') 



n ^-^ ^-^ n™ fc!(m — fc)! 



i — k\ 



(72) 



k~m — n — x' — y' + «0+ 



This series converges rapidly for large n. If Jn{N) is the 
approximation defined by summing the first N terms of 
the series (f7^ it follows that 



with 



dx 



dy 



(t>{x)(t>{y) 



n — X — y + iO+ 



(69) 



(70) 



\Jn~JniN)\ < 



\n-AK + 2\ 



max ' 



The support , [0 , — 1] , of the scaling function implies 
that in this integral x-\-y ranges from to AK — 2. This 



(73) 

where (pmax is the maximum value (< 1.5 for K — i) oi 
the scaling function. For |rt| ^ AK — 2 this error can be 
made as small as machine accuracy for modest values of 
N. 

Thus for large |n| the integrals J„ can be computed 
efficiently and accurately by truncating the sum in H72|l . 
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Integrals J„ for different values of n are related by the 
scaling equation (|64l) which when used in (|7U|I gives the 
homogeneous linear scaling equations for J„: 



Jn — hihi'J2n-l-l' 



(74) 



Equation can be used to calculate J„ recursively 
using for large |m| as input. This recursion can be 
used to step up in negative n until n = — 1 and down in 
positive n until n = AK—1. This provides an efficient and 
accurate method for calculating all of the J„ for n < 
and n > - 2. 

The remaining values, < n < 4_fC — 2, correspond to 
cases where the denominator of the singular integral 170(1 
vanishes on the support of the integrand. 

The scaling relations H74|l are still satisfied for these 
values of n, giving AK — 1 equations relating the un- 
known Jo • • • JiK~2 to the known values of J„ for n < 
and n > AK — 2. Unlike the equations for these 
equations cannot be linearly independent because they 
do not specify the treatment of the singular integral. One 
more equation is needed. 

The desired equation can be derived by observing that 
the integral J„ can be expressed in terms of the autocor- 
relation function IsOj^SlJ ^{x) of the scaling function 
as 



where 



J 11. — 



*(y) 



-dy 



(t>{x - y)(j){y)dy. 



(75) 



(76) 



It follows from the properties 

/ (j){x)dx = 1 l = ^<P{x + n) (77) 

of the scaling function that the autocorrelation function 
satisfies 

/ ^{x)dx = 1 1 = ^ $(a; + n) (78) 

and has support on [0, The autocorrelation func- 

tion is plotted in Figure 3. 

Using (|75|l in (f7B)) gives the additional linear constraint 
on the integrals J„: 



-ITT 



dx 




FIG. 3; Daubechies' K = autocorrelation function. 



which holds for any m. Replacing —in on the left side of 
equation 179|) by zero gives the principal value; by -l-iO"'" 
gives the limit on the other side of the real line. 
If m > AK — 2 equation H79|) can be expressed as 



E 



J n 



n— — m 



m+4K-3 ^4if-2 ^^^^ 



n=m+l 



-m+4:K-3 /.n+r, 

E 



$(x) 







n — x + jO+ 



dx 



dx. 



(80) 



n——m+l 

For |n| > AK — 2 the boundary integrals 



<I>(x) 



n — x + iO+ 



-dx 



<P{x) 



n — X + iO+ 



dx (81) 



and 



<i>{x) 

n — X + iO+ 



dx= I ^0+"^"^ ^^^^ 



n — X + t 



E 



E 



dx 



^{x + n) 
X - iO+ 



dx — 



(79) 



can be expanded in a convergent power series in terms of 
partial moments of the autocorrelation function 



4K-2 



<P{x) 



n — X + iO+ 



1^1 

dx = —r / ^(x)x''dx 

k=0 •'■n-m. 



(83) 
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Table 5 - J„ 



where 



rn 




J 






u 


-6.400535e 


-01 + 


iO.OOOOOOe 


J- no 


1 


-1.570288e + 00- 


i7.088321e 


m 


9 


6.615596e 


-01- 


i2.966393e 


^ uu 


Q 
O 


1.719674e + 00 + 


i6.560028e 


m 

— Ui 


4 


6.721642e 


-02- 


il.554882e 


m 

— Ui 


5 


3.595012e 


- 01 + j3.858342e 


-02 


6 


2.261977e 


-01- 


i5.023224e 


-03 


7 


1.853414e 


-01- 


i4.383938e 


-04 


8 


1.569998e 


-01- 


i3.586652e 


-06 


9 


1.357112e 


-01- 


i4.424016e 


-09 


10 


1.195044e 


-01 + 


iO.OOOOOOe 


+ 00 



and 



Jo 



$(a;) , 1^1 

-dx = — } —r 

n — X + iO' n ^-^ rv^ 

k=o 



^{x)x''dx 



(84) 

The error after truncating the series in H83f) or (|84|) after 
iV terms is bounded by 



4Jsr-2\^+^ 4:K-2 



n-AK + 2 



(85) 



where $maa: < 1 is the maximum value of the autocorre- 
lation function. It is easy to compute these quantities to 
machine accuracy. The partial moments of the autocor- 
relation function in (|83|l and (|84|) . which are needed as 
input to H80|l can be computed analytically. This calcu- 
lation is discussed in the appendix. 

Equations (|74|l and l|8()(l can be solved for the J„ for 
< n < AK — 2 in terms of left side of l|5n|l . the partial 
moments of the autocorrelation function, and the inte- 
grals for |n| > — 2. The results of these calculations 
are the nine complex numbers in Table 5. This gives val- 
ues of J„i for all TO. The solution for principal value is 
given by the real values in table 5, while the conjugate of 
the values in table gives the singular integral approaching 
the real line from the other side. 

The quantity that appears in the integral equation H49() 
is Jmn{k). In order to evaluate this quantity first note 
that it is symmetric in m and n so we can assume m > n. 
For n > we can express Jmn(k) in terms of J„; 

Jmn(k) = Jk—m — n- (^6) 

When either to or n is less than —2K + 2 then 

Jrnn{k) = (87) 

The non-trivial values of Jmn{k) correspond to to non- 
negative and -2K + 2 < n < -1, and both -2K + 2 < 
n,m < —I. This still includes an infinite number of 
integrals because k can take on any value. 

We discuss the treatment of to non-negative and 
—2K + 2 < TO < — 1 separately. When to is non-negative 
the integral becomes 



oo poo 



Um):^l dxl dy " ) (89) 
Jo J-oo m-x-y + iO+ 

Because in this case n is within 2K — 2 of zero, for large 
\m\ this can be computed in terms of moments (|23|l and 
partial moments (|24|) of the scaling function using the 
series method: 



; ^ ' k=o 

The error made by keeping N terms in the I sum is 
bounded by 



/'2K - 1 

\ m — n 



N+l 



{2K-lf 



rn — n — 2K + 1 



(91) 



which can be made as small as desired by choosing a large 
enough to. 

Using (jnH) in gives 



2K-1 2K-1 

E 

1=0 l'=0 



hihi,J2n+i{2m~l') (92) 



Jmn{k) = Jn{k - to) 



(88) 



relations among these integrals for different values of m 
and n. Equation (|92|l can be used to recursively step 
down from large values of |to| to to = —1 from below 
and TO = 2K from above. Some terms in this recursion 
will be complex because they involve the integrals J„ for 
1 < n < 4X - 3. 

The values of Jn (to) that cannot be computed directly 
from the moments or by using the recursion correspond 
-2K + 2 < n < -1 and < TO < 2ii: - 1. These can be 
computed by treating H92|l as a system of linear equations 
for the unknown J„(to)s. This works because the terms 
in (|92|l include some of the previously computed integrals. 
The results of this calculation are shown in Table 6. 

What remains are the Jmn{k) when both to and n 
fall between —2K + 2 and —1. In this case when |fc| is 
large, Jmn{k) can be expressed in the form of a conver- 
gent power series in terms of moments and partial mo- 
ments of the scaling function. The scaling equations can 
be used to step up or down in k until k is between and 
2K — 2; in addition they can also be used to solve for 
cases when 

m,n,ke [-2X+2,-l]x[-2i^+2,-l]x[0,2i4:-l]. (93) 

In a normal application the value k represents the on- 
shell energy. It will be large if there are a lot of ba- 
sis functions with support on either side of the on-shell 
point. While Jmn{k) can be calculated at the points 
from the known values using the scaling equation H92|l . 
these points do not arise in most applications. 

This completes the computation of the singular inte- 
grals that appear in equation H49|l . 
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Table 6 - J,„(n) 



in 


n 




J 




A 

-4 


U 


— o.z /Uoo4e 


— Uo H~ zu.uuuuuue 


nn 

— UU 


-4 


i 


i.ozo4yze 


— U4 — zi.oDUzyoe 


nQ 
— Uo 


A 

-4 


z 


4. izyiooe 


— U4 + ^o.yuDozie 


n/1 
— U4 


A 

-4 


Q 


— y.4(5 M / zC 


— Uo — ^y.oyziDie 


np; 
— Uo 


-4 


A 

4 


Q 1 7'Q/l /I Q„ 

o. i ( o44oe 


— UD — Zi.oUUo4ie 


— UD 


-4 







— u / — ZD. i i ( y4ze 


no 
— uy 


Q 

-o 


n 
U 


— 4.0zoooye 


— Uo H~ zu.uuuuuue 


nn 

— UU 


Q 

-o 


1 


— z.DzUOzoe 


— Uo — Z0.D4UoDie 


no 

— Uz 


Q 

-o 


n 
Z 


z.oDyzyie 


— Uz + zi.oDooooe 


no 
— Uz 


-o 


O 


— z.ou ( iUoe 


— U4 — tz.OUO ( o4e 


nQ 
— uo 


Q 


4 


— i. ioyzi4e 


— Uo — zo.ozoyiye 


n/1 
— U4 


Q 





O.4UOO0 / e 


— Uo H~ z4.y4oiioe 


— UD 




n 
U 


0. / lyozoe 


— Uz + zu.uuuuuue 


nn 
— UU 


n 
-Z 


1 


nnc 1 77„ 

o.uUDi / / e 


— Uz -j- zo.Dz4yzze 


ni 
— Ul 





o 
z 


-2.606313e 


~ 01 - i7.778517e 


no 

— Uz 


o 
-z 


Q 
O 


3.145142e 


- 02 + i2.488590e 


no 

— Uz 


-z 


A 

4 


-1.188904e 


- 02 - i5.439550e 


nQ 
— Uo 


o 
-z 





-9.599595e 


- 03 - i4.332178e 


— U4 


-1 





-2.730099e 


- 01 + iO.OOOOOOe 


-00 


-1 


1 


-4.302964e 


- 01 - il.483101e 


+ 00 


-1 


2 


1.269047e + 00 + i2.936519e 


-01 


-1 


3 


-2.565108e 


- 01 - i9.893945e 


-02 


-1 


4 


1.173431e 


- 01 + i4.010587e 


-02 


-1 


5 


4.476128e 


- 02 - i5.054775e 


-03 



V. CONCLUSION 

In this paper we introduced a method for applying 
wavelet numerical analysis to solve the relativistic three- 
body problem. The method starts by making variable 
changes in the relativistic Faddeev-Lovelace equations so 
the moving scattering singularity has simple scaling prop- 
erties. 

The next step is to project the equation in the trans- 
formed variables on a finite resolution subspace of the 
three-body Hilbert space. The matrix representation of 
the integral equation in this approximation space is easily 
computed by evaluating the driving terms and smooth 
part of the kernel at the one-point quadrature points. 
Additional integrals involving the singular part of the 
kernel over the basis functions are needed to compute the 
full kernel. These integrals can be calculated either ex- 
actly or with precisely controlled errors using the scaling 
equation and normalization condition Q. Method 
for computing all of the required integrals are discussed 
in detail in section IV and the appendix. Explicit values 
of most of the needed integrals are computed and appear 
in Tables 3-6. The physics input is in the driving term 
and smooth part of the kernel. The integrals over the sin- 
gular part of the kernel are independent of the dynamics. 
They can be computed from the scaling coefficients by 
solving a small system of linear equations. 

The resulting system of equations in the high resolu- 
tion basis, while easy to compute, is large. The wavelet 
transform is used to transform matrix elements in the 
high-resolution scaling basis to matrix elements in a basis 
consisting of low resolution scaling functions and wavelet 



basis functions with resolutions that fall between the 
high-resolution and low resolution basis. The transfor- 
mation to this new basis take 0{N') steps, which is faster 
than a fast Fourier transform. In the new basis the kernel 
can naturally be expressed as the sum of a sparse matrix 
and a small matrix. The key approximation is to replace 
the small part of the kernel by zero. The size of the error 
made in this approximation can be controlled by chang- 
ing the threshold size for discarding matrix elements. 

The sparse matrix can be solved using sparse ma- 
trix techniques. In reference ^ this was done by first 
conditioning the matrix and then using the complex bi- 
conjugate gradient method. The resulting solution can 
be transformed back to the scaling basis using the in- 
verse wavelet transform. The approximation is improved 
if the resulting solution is substituted back in the original 
equation (2^. The step has the added benefit that the 
basis functions never have to be calculated. 

The key result that was needed to calculate integrals 
associated with the mo ving singularities is the observa- 
tion by Beylkin [2^ [s^lllj and collaborators that inte- 
grals of scaling functions over moving singularities can be 
expressed as integrals of the autocorrelation function of 
the scaling function over a fixed singularity. This leads 
to a practical and stable method for computing the inte- 
grals. The methods do not require subtractions or careful 
choices of quadrature points; for the Daubechies' K = Z 
basis they are reduced to solving a system of eleven linear 
equations. The required properties of the autocorrelation 
function are derived in the appendix. 

The research in references [l| and demonstrated 
that the wavelet method led to sparse matrix approxima- 
tions resulting negligible errors. The structure of the ker- 
nel in the relativistic three-body case indicates that the 
wavelet method will lead to accurate sparse matrix ap- 
proximations to the relativistic Faddeev-Lovelace equa- 
tions. 

The increase in eflaciency in this method is due to the 
saving in computational effort in going from solving a 
large dense set of linear equations to an approximately 
equivalent equations with a sparse matrix. The wavelet 
method will lead to a significant savings in computational 
effort for a large system. 

Our conclusion is that wavelet numerical analysis can 
be used to accurately approximate the relativistic Fad- 
deev Lovelace equations by a linear system of equations 
with a sparse kernel matrix. 
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Table 7: Daubechies' K — 3 Autocorrelation Scaling 
Coefficients 





7.825529e 


- 02 


ai 


3.796160e 


- 01 


0-2 


6.767361e 


- 01 


0.3 


4.612557e 


- 01 


04 


-4.471656e 


-02 


as 


-1.687321e 


-01 


06 


-2.481571e 


-03 


07 


3.922363e 


-02 


08 


-1.563883e 


-03 


Og 


-4.256471e 


-03 


OlO 


8.774429e 


-04 



APPENDIX A 

The autocorrelation function of the scahng function is 
defined by 

/oo 
cl>{x - y)cb{y)dy. (Al) 
-oo 

Because of the support of the scahng function is [0, 2K ~ 
1], the autocorrelation function has support [0,4i4r — 2]. 

Using the scaling equation Ht)4|) for the scaling func- 
tion in the definition ljAl|l of the autocorrelation func- 
tion leads to the scaling equation for the autocorrelation 
function: 

$(x) = ^^/i,/ip$(2a;-Z'-0- (A2) 
If we define 

ai = ^ hi-i'^i' 0<l<4K-2 (A3) 

equation IIA2p can be put in the same form as 

iK-2 

D<i>{x) = aiT^<^{x). (A4) 

1=0 

The scaling coefficients ak for the Daubechies' K — 3 
autocorrelation function are given in Table 7. 



The normalization condition of the scaling function 
can be used in the definition IjAip of the autocorrelation 
function to derive the normalization condition: 

J <^>{x)dx = 1. (A5) 

The scaling equation ljA2|l and the normalization con- 
dition (|A5() can be used calculate moments and partial 



Table 8: Daubechies' K — 3 Autocorrelation 
Moments 

(x-")s l.OOOOOOe + 00 

(s^)* 1.634802e + 00 

(x^)* 2.672579e + 00 

(x^)* 4.167773e -hOO 

(x*)* 5.825913e + 00 

(x^)* 6.817542e + 00 

(x'^)s 8.8079176 4-00 

{x^)<i, 4.055470e + 01 

{x*)s 2.899550e + 02 

{x'')<i, 1.695851e + 03 

(x^°).i, 8.321402e -I- 03 



moments of the autocorrelation function. To calculate 
the moments of the autocorrelation function use 

(x*^)* = / <^>{x)x''dx, (A6) 



Moving the n — k term to the left side of the equation 
gives recursion relation 

(A8) 

The recursion is started using the the normalization con- 
dition (|A5p . The lowest moments are tabulated in Table 
8. 

The partial moments of the autocorrelation function 
satisfy the scaling equation 

(a;'')<i,,[,„,oo] := / <i>{x)x''dx = 



I n=0 

When m > 4K — 2 these partial moments become ordi- 
nary moments while when m < they vanish. This gives 
us a linear system for partial moments in terms of the 
full moments and lower partial moments. These equa- 
tions can be solved recursively. Partial moments corre- 
sponding to more general intervals can be computed by 
subtraction: 

(x'')<i>,[m,n] = (a^'')$,[m,oo] " {x'') <i> ,[n,oo] ■ (AlO) 
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